Optosurf simulation
1 Gaussian function definition
Code
# Create a function to plot the equation
def plot_equation(mu, sigma, n, number_points, degrees, plot, title="Super-Gaussian", width = 700, height = 550):
"""
Plot the Super-Gaussian equation using Bokeh
Parameters
----------
mu (float): Mean value for the equation
sigma (float): Standard deviation value for the equation
n (float): Order value for the equation
number_points (int): number of points to calculate the function
Returns
-------
plots (bokeh plot): Plot of the Super-Gaussian equation
x(np): linspace for the gaussian plot
y(np): gaussian values
"""
# 1. Define linear degrees vector and calculate Super-Gaussian
ticker = SingleIntervalTicker(interval=2.5, num_minor_ticks=10)
xaxis = LinearAxis(ticker = ticker)
x = np.linspace(degrees[0], degrees[1], number_points)
y = np.exp(-abs(((x-mu)/sigma))**n)
# 2. Plot
if plot:
TOOLTIPS = [("index", "$index"),("(x,y)", "($x, $y)")]
p = figure(title=title, x_axis_label='x', y_axis_label='y', tooltips = TOOLTIPS,
width = width, height = height)
p.line(x, y, line_width=4, alpha = 0.5, line_color = "#C5E064")
p.add_layout(Grid(dimension=0, ticker=xaxis.ticker))
p = plot_format(p, "Degrees", "Intensity", "bottom_left", "10pt", "10pt", "10pt")
return p, x, y
else:
return x, y
p, x, y = plot_equation(0, 1, 3.5, 50000, [-15, 15], True)
show(p)1.1 Window integration
The optosurf head has a 32 pixels line detectors. This is simulated by performing a window integration:
Code
# Create a function to do the window integration
def window_integration(number_windows, window_size, gap, x, y, p=None):
"""
Performs a window integration
Parameters
----------
number_windows (int): Number of integration windows
window_size (int): Number of data points in the window
x(np): linspace for the gaussian plot
y(np): gaussian values
Returns
-------
p (bokeh plot): Plot of the integration
integration_axis (np): window integration axis
integration_points (np): Integrated points
"""
integration_points = []
integration_axis = []
color_multiplier = len(bokeh.palettes.Viridis256)//number_windows
count = 0
for i in range(number_windows):
# 1. Get data in every window and integrate
a = i*window_size
b = i*window_size + window_size
x_temp = x[a:b-gap:1]
y_temp = y[a:b-gap:1]
integration = np.trapz(y_temp, x_temp, dx = x[1] - x[0])
integration_points.append(integration)
axis = x_temp[len(x_temp)//2]
integration_axis.append(axis)
# 2. Draw a rectangle of the integration window
if p is not None:
left_edge = x_temp[0]
right_edge = x_temp[-1]
p.rect(x=(left_edge + right_edge)/2, y=0.18, width=right_edge-left_edge, height=0.3, fill_alpha=0.001, fill_color='#C5E0B4', color='#C5E0B4')
p.rect(x=(right_edge + x[b-1])/2, y=0.18, width=x[b-1]-right_edge, height=0.3, fill_alpha=0.005, fill_color='#F16C08', color = '#F16C08')
p.circle(x_temp[::15], y_temp[::15], size = 4, alpha = 1)
count += 1
if p is not None:
p.circle(integration_axis, integration_points, size = 7, color = '#FAA0A0')
p.line(integration_axis, integration_points, line_width = 4, color = '#FAA0A0', alpha = 0.8)
return p, integration_axis, integration_points
p, int_axis, int_points = window_integration(32, 1562, 100, x, y, p)
show(p)The optosurf head has a 32 pixels line detectors. This is simulated by performing a window integration:
2 Histogram
The next step is to make a histogram reconstruction to calculate the standard deviation
Code
# Create a function to do histogram reconstruction
def histogram_reconstruction(int_points, hist_bool):
"""
Constructs a histrogram
Parameters
----------
int_points(np): Points calculated from the window integration
Returns
-------
hist_plot (bokeh plot): Plot of the histogram
std_dev(float): Histogram's standard deviation
"""
# a. Histogram reconstruction
normalized_y = np.multiply(int_points, 10000)
hist_2d = np.array([])
for i, int_point in enumerate(normalized_y):
round_int_point = round(float(int_point))
hist_2d = np.concatenate((hist_2d, np.array([float(i)]*round_int_point)))
# b. Calculate standard deviation
stddev = np.std(hist_2d)
# c. Plot histogram
if hist_bool:
hist_plot = alt.Chart(pd.DataFrame({'x': hist_2d})).mark_bar(opacity=0.7, color='#2D908C').encode(
alt.X('x:Q', bin=alt.Bin(maxbins=20)),
y='count()',
tooltip=['x','count()'])
title = f'Standard deviation: {stddev}'
hist_plot = hist_plot.properties(title=title, background='#282B30')
return hist_plot, np.std(stddev)
else:
return stddev
alt.data_transformers.disable_max_rows()
alt.renderers.set_embed_options(theme='dark')
hist_plot, std_dev = histogram_reconstruction(int_points, True)
hist_plot3 Standard deviation matrix
Code
mu_np = np.linspace(-5, 5, 4)
std_np = np.linspace(1.0, 5, 7)
X, Y = np.meshgrid(mu_np, std_np)
std_grid = np.empty_like(X)
plots_gaussian = []
n = 3.5
number_points = 10000
degrees = [-15, 15]
number_windows = 32
window_size = number_points//number_windows
gap = 100
for i in range(len(mu_np)):
for j in range(len(std_np)):
# Generate x and y Gaussian data points
title = f"mu: {mu_np[i]:.3f}, std: {std_np[j]:.3f}"
# print(title)
plot, x, y = plot_equation(mu_np[i], std_np[j], n, number_points, degrees, True, title, 200, 200)
plot.title.text_font_size = "10pt"
plot, int_axis, int_points = window_integration(number_windows, window_size, gap, x, y, plot)
plots_gaussian.append(plot)
grid_gaussian = gridplot(children = plots_gaussian, ncols = len(std_np), merge_tools=False)
show(grid_gaussian)